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ABSTRACT 

o 

Context. Classical novae are explosive phenomena that take place in stellar binary systems. They are powered by mass transfer from 
a low-mass, main sequence star onto a white dwarf. The material piles up under degenerate conditions and a thermonuclear runaway 
ensues. The energy released by the suite of nuclear processes operating at the envelope heats the material up to peak temperatures of 
~ (1 -4) x 10 8 K. During these events, about 10~ 4 - 10~ 5 M©, enriched in CNO and other intermediate-mass elements, are ejected into 
the interstellar medium. To account for the gross observational properties of classical novae (in particular, a metallicity enhancement 
in the ejecta above solar values), numerical models assume mixing between the (solar-like) material transferred from the companion 
and the outermost layers (CO- or ONe-rich) of the underlying white dwarf. 

Aims. The nature of the mixing mechanism that operates at the core-envelope interface has puzzled stellar modelers for about 40 years. 
Here we investigate the role of Kelvin-Helmholtz instabilities as a natural mechanism for self-enrichment of the accreted envelope 
with core material. 

Methods. The feasibility of this mechanism is studied by means of the multidimensional code FLASH. Here, we present a series 
of 9 numerical simulations perfomed in two dimensions aimed at testing the possible influence of the initial perturbation (duration, 
strength, location, and size), the resolution adopted, or the size of the computational domain on the results. 

Results. We show that results do not depend substantially on the specific choice of these parameters, demonstrating that Kelvin- 
Helmholtz instabilities can naturally lead to self-enrichment of the accreted envelope with core material, at levels that agree with 
observations. 

Key words. (Stars:) novae, cataclysmic variables — nuclear reactions, nucleosynthesis, abundances — convection — hydrodynamics 
— instabilities — turbulence 
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| 1. Introduction cal simulation of such outbursts. To match the energetics, peak 

. luminosities, and the abundance pattern, models of these explo- 

} Classical novae are cataclysmic stellar events. Their ther- sio ns require mixing of the material accreted from the low-mass 
f x | monuclear origin, theorized by Schatzmann (1949, 1951) and steUar companion with the outer layers of the underlying white 
Cameron (1959) - see also Gurevitch & Lebedinsky (1957) dwarf In fact> because of the moderate temp eratures achieved 
© . and references therein - has been established through multi- during the TNR a very Hmited production of elements beyond 
t^h wavelength observations and numerical simulations pioneered those from the C NO-cycle is expected (Starrfield et al. 1998, 
> by Sparks (1969), who performed the first 1-D, hydrodynamic 2009; ]qs6 & Hernanz 199g; Kovetz & PriaMk 199?; Yaron et 
nova simulation. These efforts helped to establish a basic pic- al 2005)? and the ific chemical abu ndances derived from ob- 
i ture, usually referred to as the thermonuclear runaway model servations (with a suite f elements ranging from H to Ca) cannot 
rrt (TNR), which has been successful m reproducing the gross ob- be explained by thermonuclear processing of solar-like material, 
servational properties of novae, namely the peak luminosities Mixing at the core _ envelope interface represents a likely mecha- 
achieved, the abundance pattern, and the overall duration of the n ism 
event; see Starrfield et al. (2008), Jose & Shore (2008), Jose & 

Hernanz (2007) for recent reviews. The details of the mixing episodes by which the envelope 

Many details of the dynamics of nova explosions remain is enriched in metals have challenged theoreticians for nearly 
to be explored. In particular, there are many observed cases of 40 years. Several mechanisms have been proposed, including 
nonspherical ejecta, inferred from line profiles during the early diffusion-induced mixing (Prialnik & Kovetz 1984; Kovetz & 
stages of the outburst and from imaging of the resolved ejecta, Prialnik 1985; Iben et al. 1991, 1992; Fujimoto & Iben 1992), 
including multiple shells, emission knots, and chemical inhomo- shear mixing at the disk-envelope interface (Durisen 1977; 
geneities. Although the broad phenomenology of the outburst Kippenhahn & Thomas 1978; MacDonald 1983; Livio & Truran 
can be captured by 1-D calculations, it is increasingly clear that 1987; Kutter & Sparks 1987; Sparks & Kutter 1987), convec- 
the full description requires a multidimensional hydrodynami- tive overshoot-induced flame propagation (Woosley 1986), and 

mixing by gravity wave breaking on the white dwarf surface 

Send offprint requests to: J. Jose (Rosner et al. 2001; Alexakis et al. 2004). The multidimen- 
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sional nature of mixing has been addressed by Glasner & Livne 
(1995) and Glasner et al. (1997, 2005, 2007) with 2-D sim- 
ulations of CO-novae performed with VULCAN, an arbitrar- 
ily Lagrangian Eulerian (ALE) hydrocode capable of handling 
both explicit and implicit steps. They report an effective mixing 
triggered by Kelvin-Helmholtz instabilities that produced metal- 
licity enhancements to levels in agreement with observations. 
Similar studies (using the same initial model as Glasner et al. 
1997) were conducted by Kercek et al. (1998, 1999) in 2-D and 
3-D, respectively. Their results, computed with the Eulerian code 
PROMETHEUS, displayed mild TNRs with lower peak temper- 
atures and velocities than Glasner et al. (1997) and insufficient 
mixing. While Glasner et al. (1997) argue that substantial mixing 
can naturally occur close to peak temperature, when the envelope 
becomes fully convective and drives a powerful TNR, Kercek et 
al. (1998) conclude instead that mixing must take place much 
earlier: if it occurs around peak temperature, it leads to mild ex- 
plosions or to events that do not resemble a nova. 

The differences between these studies have been carefully 
analyzed by Glasner et al. (2005), who conclude that the early 
stages of the explosion, before TNR ignition when the evolu- 
tion is quasi- static, are extremely sensitive to the adopted outer 
boundary conditions. They show that Lagrangian simulations, 
in which the envelope is allowed to expand and mass is con- 
served, lead to consistent explosions. In contrast, in Eulerian 
schemes with a "free outflow" outer boundary condition, the 
choice adopted in Kercek et al. (1998), the outburst can be arti- 
ficially quenched. The scenario was revisited by Casanova et al. 
(2010), who show that simulations with an Eulerian scheme — 
the FLASH code — and a proper choice of the outer boundary 
conditions can produce deep-mixing of the solar-like accreted 
envelopes with core material. The puzzling results reported by 
Kercek et al. (1998) stress the need for a systematic evaluation 
of the effect that different choices of model parameters (e.g. the 
intensity and location of the initial temperature perturbation, res- 
olution, or size of the computational domain) may have on the 
results. To this end, we performed a series of 9 numerical simu- 
lations in 2-D aimed at testing the influence of these parameters 
on the level of metal enhancement of the envelope. Here we re- 
port the results of these simulations. 

Our paper is organized as follows. In Sect. 2 we explain our 
input physics and initial conditions. Then Sect. 3 is devoted to 
studying the mixing at the core-envelope interface for our fidu- 
cial model. In Sect. 4 the effects of the size of the initial pertur- 
bation are analyzed, while in Sect. 5 we discuss the effects of 
the size of the computational domain. In Sect. 6 we quantify the 
influence of the grid resolution. Finally, in Sect. 7 we discuss the 
significance of our results and draw our conclusions. 

2. Input physics and initial conditions 

The simulations reported here were performed with FLASH, 
a parallelized explicit Eulerian code, based on the piecewise 
parabolic interpolation of physical quantities for solving the hy- 
drodynamical equations, and with adaptive mesh refinement (see 
Fryxell et al. 2000). As in Casanova et al. (2010), we used the 
same initial model as Glasner et al. (1997) and Kercek et al. 
(1998): a 1 Mq CO white dwarf that accretes solar composition 
matter (Z = 0.02) at a rate of 5 x 10" 9 M Q yr" 1 . The model was 
evolved spherically (1-D) and mapped onto a 2-D cartesian grid, 
when the temperature at the base of the envelope reached « 10 8 
K. It initially comprised 112 radial layers — including the outer- 
most part of the CO core — and 512 horizontal layers. The mass 
of the accreted envelope was about 2 x 10" 5 M©. Nuclear energy 
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Fig. 1. Snapshot of the development of early instabilities, which 
later spawn Kelvin-Helmholtz instabilities, shown in terms of 
the 12 C mass fraction (in logarithmic scale) for model A, 158 s 
from the start of the simulation when the core-envelope interface 
temperature is 7b aS e ~ 1.36 x 10 8 K. 

generation is handled through a network of 13 species (*H, 4 He, 

12,13 Cj 13,14,15 Nj 14,15,16,17 0j and 17,18p )? and connected through 

18 nuclear reactions. We adopted the conductive and radiative 
opacities from Timmes (2000) and an equation of state based on 
Timmes & Swesty (2000). Periodic boundary conditions were 
imposed on both sides of the computational domain with ver- 
tical hydrostatic equilibrium with an outflow constraint at the 
top and a reflecting constraint at the bottom on the velocity (see 
Zingale et al. 2002). A summary of the main characteristics of 
the 9 models computed in this work is given in Table 1 , where H 
is the distance from the perturbation to the initial core-envelope 
interface, R x and ST, and St are the size, strength, and du- 
ration of the temperature perturbation, and Z the mass-averaged 
metallicity of the envelope at the end of the calculations. 

3. 2-D simulations of mixing at the core-envelope 
interface 

In this Section, we describe the basic features of our fiducial 
model A, as a framework for further discussion of the effect of 
the parameter choices on our results. A movie, showing the de- 
velopment of Kelvin-Helmholtz instabilities, in terms of the 12 C 
content, up to the time when the convective front hits the up- 
per computational boundary, ModelA-2D.wmv, is available at 
http://www.fen.upc.edu/users/jjose/Downloads.html The simu- 
lation was performed for the conditions of model A, as summa- 
rized in Table 1 . 

For all sequences reported in this work, the relaxation of the 
initial model to guarantee hydrostatic equilibrium, together with 
the small amount of numerical viscosity — in contrast with the 
simulations performed by Glasner et al. (1997) — requires an 
initial perturbation close to the core-envelope interface to trig- 
ger the onset of instabilities early in the calculations. The initial 
perturbation is applied to the temperature using four parame- 
ters: strength, location, size and duration. Model A assumes a 

1 The different values adopted for R x and R y in models F and G are 
used to account for the assumption of a rectangular (rather than square) 
computational domain. 



Table 1. Models computed. 
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Fig. 2. Mach number at two different moments of the simulation, t = 230 s (left panel) and 496 s (right panel), for model A. 



top-hat temperature perturbation wherever ((x - xq) /R x ) 2 + ((y - 
yo)/R y ) 2 < 1, where x and y are the space coordinates measured 
from the center of the perturbation, (jto, yo), and R x and R y indi- 
cate its spatial extent. We fixed xq = 5 x 10 7 cm in all sequences. 
The strength of the perturbation is 5% in temperature in all cases 
but one (see table 1). It is 1 km wide, applied only during the ini- 
tial timestep (that is, the temperature is fixed only during 10" 10 
s), and imposed on the core-envelope interface (yo = 5.51 x 10 8 
cm). The resolution adopted in model A is 1.56x 1.56 km, and 
the size of the computational domain is 800 x 800 km. 

The initial perturbation drives a shear flow that triggers the 
formation of instabilities (Fig. 1), about 150 s after the start of 
the simulation. As soon as material from the core is mixed into 
the envelope, small convective cells develop. At this early stage, 
the fluid has a large Reynolds number, with a characteristic eddy 
length of 50 km, fluid velocities ranging between v = 10 5 - 10 6 
cm s" 1 , and a dynamic viscosit>0 of 2 x 10 4 P The fluid ve- 
locity v remains below the speed of sound c s (that is, the Mach 
number Ma= v/c s is always less than unity, see Fig. 2), hence, 
the fluid displays a deflagration rather than a detonation — see 
Williams (1985), for a thorough analysis of the differences be- 
tween flame propagation under detonation and deflagration con- 
ditions. At t = 235 s, The characteristic eddy turnover time is 
^conv/vconv ~ 10 s. The merging of the small convective cells into 
large eddies, characteristic of 2-D simulations, with a size com- 



2 The dynamic viscosity evaluates the resistance to flow of a fluid 
under an applied force. More precisely, it is defined as the tangential 
force per unit area required to move one horizontal plane with respect 
to the other at unit velocity when maintaining a unit distance apart by 
the fluid. 



parable to the height of the envelope, reinforces the injection of 
CO-rich material into the envelope. Convection becomes more 
turbulent. At this stage (t = 450 s), the nuclear energy genera- 
tion rate exceeds 10 15 erg g" 1 s" 1 , while the characteristic con- 
vective timescale decreases to ~ 5 s. The convection front prop- 
agates progressively upwards (Fig. 3, left panel), with a velocity 
of ~ 10 km/s, and eventually reaches the top of our computa- 
tional domain. The envelope base reaches a peak temperature of 
1.64 x 10 8 K. At this time (t = 496 s), when matter starts to 
cross the outer boundary of the computational domain, we stop 
the calculations because of the Eulerian nature of the FLASH 
code. At this final stage, the mean mass-averaged metallicity 
in the envelope reaches Z ~ 0.21. It is worth noting, however, 
that the convective eddies are still pumping metal-rich matter 
through the core-envelope interface. Hence, it is likely that the 
final metallicity in the envelope will be larger. The simulation 
shows that the induced Kelvin-Helmholtz vortices can naturally 
lead to self-enrichment of the accreted envelope with core mate- 
rial to levels that agree with observations and that the expansion 
and progress of the runaway is almost spherically symmetric for 
nova conditions even for a point-like TNR ignition. 



4. Effect of the initial perturbation 

To quantify the influence of the initial perturbation on our re- 
sults, we have performed a series of 2-D hydrodynamic tests for 
a set of different durations, strengths (intensities), locations and 
sizes of the perturbation. For simplicity, a top-hat perturbation, 
centered at xo = 5 x 10 7 cm, has been adopted in all models 
reported in this work. 
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Fig. 3. Left panel: propagation of the convective front as a function of time, for models A, H, and I. Right panel: temperature profile 
versus radius at two different times, t = s (solid line; Tbase = 9.84 x 10 7 K) and t = 496 s (dashed line; 7b aS e = 1.64 x 10 8 K), for 
model A. 



The effect of the duration of the perturbation was checked by 
means of a test case (model B), identical to model A but with a 
perturbation lasting for 10 s. As shown in Table 1, the charac- 
teristic timescales for model B, such as the time required for the 
first instabilities to show up, Tkh, or the time needed by the con- 
vective front to hit the outer boundary, ty, become shorter. The 
role played by a temperature perturbation can be understood in 
terms of the energy injected into the envelope: the longer the 
duration of the perturbation, the larger the energy injected, and 
thus, the shorter the characteristic timescales of the TNR. This 
has little effect, however, on the overall metallicity enhancement 
in the envelope since a final CNO mass fraction of ~ 0.212 was 
found in model B, whereas ~ 0.224 resulted in model A. 

Both models A and B assumed temperature perturbations of 
ST ~ 5% during the initial timestep (~ 10" 10 s). To test the 
possible influence of the strength of the perturbation, a test case 
with ST ~ 0.5% (model C) has also been computed. As shown in 
Table 1 and Fig. 4, the time evolution of models A and C is very 
similar, and hence, similar final mean CNO mass fractions at the 
end of the simulations were found (with Z = 0.209 in model C). 

The effect of the location of the perturbation along the ver- 
tical axis has also been studied: whereas model A assumed a 
temperature perturbation of ~ 5%, applied at the innermost en- 
velope shell (yo = 5.51 x 10 8 cm), in model D, a similar per- 
turbation was placed ~ 5 km above the core-envelope interface 
(yo = 5.5 15 x 10 8 cm). Both models exhibit a very similar tempo- 
ral evolution, with almost identical times for the appearance of 
the first instabilities and for the time required to reach the outer 
boundary. Similar envelope mean CNO mass fractions (0.224 
and 0.235, respectively) were also found. 

Finally, the influence of the size of the perturbation has also 
been analyzed. Whereas model D was evolved with an initial 
temperature perturbation of size R x = 1 km and R y = 1 km, model 
E assumed R x = 5 km and R y = 5 km. As before, very similar 
characteristic timescales (see Table 1) and final mean CNO mass 
fractions (0.235 and 0.209, respectively) were found. 

To summarize, the specific choice of the parameters that de- 
fine the initial temperature perturbation has a negligible effect on 
metallicity enhancement of the envelope. 

5. Effect of the size of the computational domain 

The choice of the computational domain represents a compro- 
mise between computational time requirements and numerical 



accuracy. Several considerations constrain its minimum size. On 
one hand, the merger of large convective eddies often found in 
2-D simulations may be severely affected by the adoption of a 
small computational domain. On the other hand, nova outbursts 
eventually result in mass ejection. With an Eulerian code such as 
FLASH, it is not possible to track the material that flows off the 
grid, and hence, it is important to use domains that are as large 
as possible along the radial direction (while being sufficiently 
wide along the horizontal axis). Unfortunately, when the initial 
1-D model is mapped into the 2-D grid, and relaxed to guarantee 
hydrostatic equilibrium, densities quickly underflow values for 
large heights (Zingale et al. 2002). 

The specific size adopted for most of the models computed 
in this work, i.e. 800 x 800 km, is a bit smaller than those used in 
Glasner et al. (1997) — 0.l7r rad , in spherical-polar coordinates 
— and in Kercek et al. (1998) — 1800 x 1 100 km, in a cartesian, 
plane-parallel geometry. In this section, we analyze possible de- 
pendences of the results on the adopted size of the computational 
domain. To this end, two additional simulations were performed. 
In the first one (model F), a wider computational domain has 
been adopted (i.e., 1600 x 800 km). In the second (model G), 
aimed at testing the influence of the vertical (radial) length, a 
domain of 800 x 1000 km has been used (where the choice of 
1000 km results from numerical restrictions that limit the verti- 
cal extent of our computational domain). 

As shown in Table 1 , the horizontal width (model F) has no 
noticeable effect on the timescales of the simulations, either for 
the time required for the onset of the first instabilities or for the 
time required for the convective front to reach the outer bound- 
ary. The mass-averaged CNO abundance in the envelope reached 
~ 0.206 at the end of this simulation, close to the value found 
for model A. These results confirm that 800 km is an appropri- 
ate choice for the width of the computational domain, stressing 
that above a threshold value the course of the TNR is insensi- 
tive to the adopted width, in agreement with the sensitivity study 
performed by Glasner et al. (2007). 

The specific length adopted along the vertical direction (see 
model G), while unimportant for the time of appearance of the 
instabilities (around 155 s after the start of the simulation, as in 
model A), affects the time required to reach the outer boundary, 
located 200 km above the value adopted for model A. Moreover, 
the larger extension of the computational domain along the radial 
(vertical) direction allows the convective eddies to pump addi- 
tional metal-rich core material into the envelope compared with 
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Fig. 4. Left panel: time evolution of the nuclear energy generation rate (in erg.s l ) for the 9 models computed in this work. Right 
panel: final CNO mass fraction versus radius. 



all the simulations reported previously in this paper. Indeed, the 
mean, mass-averaged metallicity in model G achieves the largest 
value of all the simulations reported, ~ 0.291. This result sug- 
gests that the likely mean mass-averaged metallicity driven by 
Kelvin-Helmholtz instabilities should be Z « 0.3. In summary, 
we conclude that the size of the computational domain, above a 
certain threshold value, has little influence on the physical quan- 
tities that are more directly related with the mixing process at the 
core-envelope interface. 



6. Effect of the grid resolution 

All simulations discussed so far (e.g., models A to G) were per- 
formed with a resolution of 1.56 x 1.56 km, a value similar to 
the minimum resolutions adopted in Glasner et al. (2007) which 
is roughly ~ 1.4 x 1.4 km, and in Kercek et al. (1998), 1 x2 km. 
To quantitatively assess the possible effect of the resolution, two 
additional test cases were computed with exactly the same input 
parameters as in model A but with two different resolutions: lxl 
km (model H) and 0.39 x 0.39 km (model ifl 

As shown in Table 1, the increase in resolution produces a 
delay in the time required for the first instabilities to develop, 
£kh- This seems to be a numerical artifact. In models with a 
coarser resolution, the larger size of the blocks artificially gen- 
erates a larger numerical diffusion compared to models with a 
finer resolution (a similar resolution dependence is clearly seen 
as well in the Kercek et al. 1998 simulations). Actually, the ratio 
of differences in the initial build up times (i.e., (model I-model 
A)/(model H-model A)) scales approximately as the zone size 
dimensions to the power of two. This is a purely numerical per- 
turbation that forces the development of instabilities. To test this 
hypothesis, we computed an additional test case (not included 
in Table 1), identical to model A but without any initial pertur- 
bation. The onset of the instabilities in such an extremely low 
numerical diffusion regime is substantially delayed. The simula- 
tions reported by Glasner et al. (1997) also show the early ap- 
pearance of instabilities in a model with substantial numerical 
noise: within a very short time (about 10 s), the numerical noise 



3 For comparison, whereas a maximum number of 5300 blocks are 
administered in model A, the number of blocks increases up to 83000 
in model I. The total CPU time spent in both simulations, using 256 
processors of the MareNostrum supercomputer, has been 3 and 110 khr, 
respectively. 



( round-off) seeds an intense convectiveflow in the envelope with- 
out any artificial perturbations. 

A similar behavior is also found for the time required for the 
convective front to reach the outer boundary, ty, and for the his- 
tory of the nuclear energy generation rate (Fig. 4). As expected, 
filamentary structures and convective cells are better resolved in 
the finer resolution model I, compared to those computed with 
somewhat coarser grids (models A and H; see Fig. 5). These mi- 
nor differences do not, however, show significant variations in 
the final, mean CNO abundances achieved in the envelope: while 
Z ~ 0.224 in model A, models H and I yield 0.201 and 0.205, by 
mass, respectively. Similar agreement is found in the peak tem- 
peratures achieved and in the overall nuclear energy generation 
rates (Fig. 4). 

Thus, the adopted resolution has not a critical effect for the 
mixing models presented in this work. The variation in the final 
mean CNO abundance in the envelope, under the range of reso- 
lutions adopted, is only about 12% (when comparing results for 
models A, H, and I), a quite reasonable value. 



7. Discussion and Conclusions 

In this paper we have reported results for a series of nine 2-D 
numerical simulations that test the influence of the initial pertur- 
bation (duration, strength, location, and size), the resolution of 
the grid, and the size of the computational domain on the results. 
We have shown that mixing at the core-envelope interface pro- 
ceeds almost independently of the specific choice of such initial 
parameters, above threshold values. 

The study confirms that the metallicity enhancement inferred 
from observations of the ejecta of classical novae can be ex- 
plained by Kelvin-Helmholtz instabilities, powered by an ef- 
fective mesoscopic shearing resulting from the initial buoyancy. 
Fresh core material is efficiently transported from the outermost 
layers of the white dwarf core and mixed with the approximately 
solar composition material of the accreted envelope. As soon as 
12 C and 16 are dredged up, convection sets in and small con- 
vective cells appear, accompanied by an increased nuclear en- 
ergy generation rate. The size of these convective cells increases 
in time. Eventually, these cells merge into large convective ed- 
dies with a size comparable to the envelope height. The range 
of mean mass-averaged envelope metallicities obtained in our 
simulations at the time when the convective front hits the outer 
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Fig. 5. Snapshots of the l H (upper panels) and 12 C (middle panels) mass fractions at t ~ 395 s (model A; left panels), and 688 s 
(model I; right panels). Lower panels: the number of blocks administered, at this stage, is 3184 for model A, and 43800 for model 
I. In both simulations, FLASH divides each block in 8 cells. Structures such as vortexs are better resolved in the finer resolution 
model I. 



boundary, 0.21 - 0.29, matches the values obtained for classical 
novae hosting CO white dwarfs. 

It is, however, worth noting that the convective pattern is 
actually produced by the adopted geometry (e.g., 2-D), forcing 
the fluid motion to behave very differently than 3-D convection 
(Shore 2007; Meakin & Arnett 2007). Nevertheless, the levels of 
metallicity enhancement found in our 2-D simulations will likely 
remain unaffected by the limitations imposed by the 2-D geom- 



etry (D. Arnett, private communication). Fully 3-D simulations 
aimed at testing this hypothesis are currently underway. 
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